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Abstract. We present the first automated fitting method for the quantitative spectroscopy of O- and early B -type stars with 
stellar winds. The method combines the non- LTE stellar atmospher e code fastwmd from lPuls et alJ (2005) with the genetic 
algorithm based optimizing rou tine pikaia fr om Charbonneau 1 1995), allowing for a homogeneous analysis of upcoming large 
samples of early-type stars ("e.g. lEvans et all2 0051. In this first implementation we use continuum normalized optical hydrogen 
and helium lines to determine photospheric and wind parameters. We have assigned weights to these lines accounting for line 
blends with species not taken into account, lacking physics, and/or possible or potential problems in the model atmosphere 
code. We find the method to be robust, fast, and accurate. Using our method we analysed seven O-type stars in the young 
cluster Cyg OB2 and five other Galactic stars with high rotational velocities and/or low mass loss rates (including 10 Lac, 
£ Oph, and t Sco) that have been studied in detail with a previous version of fastwind. The fits are found to have a quality that 
is comparable or even better than produced by the classical "by eye" method. We define errorbars on the model parameters 
based on the maximum variations of these parameters in the models that cluster around the global optimum. Using this concept, 
for the investigated dataset we are able to recover mass-loss rates down to ~ 6 x 10~ 8 MQyr~' to within an error of a factor 
of two, ignoring possible systematic errors due to uncertainties in the continuum normalization. Comparison of our derived 
spectroscopic masses with those derived from stellar evolutionary models are in very good agreement, i.e. based on the limited 
sample that we have studied we do not find indications for a mass discrepancy. For three stars we find significantly higher 
surface gravities than previously reported. We identify this to be due to differences in the weighting of Balmer line wings 
between our automated method and "by eye" fitting and/or an improved multidimensional optimization of the parameters. The 
empirical modified wind momentum relatio n constru cted on the basis of the stars analysed here agrees to within the error bars 
with the theoretical relation predicted bv lVinketalJ r2000). including those cases for which the winds are weak (i.e. less than 
a few times 10~ 7 Moyr -1 ). 

Key words, methods: data analysis - line: profiles - stars: atmospheres - stars: early-type - stars: fundamental parameters - 
stars: mass loss 



1. Introduction 

Until about a decade ago detailed analysis of the photospheric 
and wind properties of O-type stars was limited to about 40 to 
50 stars divid ed over the Galaxy a nd the M agellan ic Clouds 
(see e .g. IPuls et al1ll996l see also iReoolust. Puis. & Herrerol 
120041) . The reason that at that time only such a limited number 
of objects had been investigated is related in part to the fact that 
considerable effort was directed towards improving the physics 
of the non-local thermodynamic equilibrium (non-LTE) model 
atmospheres used to analyse massive stars. Notable devel- 
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opm ents have been the im provements in the atomic models 
(e.g. [Becker & Buflerlll992l). shock treatment JPauldrach et alJ 
1200 lb . clumping failherll99ll:lHillier & Millerll999l). and the 
implementation of line blanketing (e.g.lHub env & L anzll 19951 
Iffillier & MiUerlll998tlPauidrach et alJbOOlh . To study the ef- 
fects of these new physics a core sample of "standard" O-type 
stars has been repeatedly re-analysed. A second reason, that 
is at least as important, is the complex, and time and CPU 
intensive nature of these quantitative spectroscopic analyses. 
Typically, at least a six dimensional parameter space has to 
be probed, i.e. effective temperature, surface gravity, helium 
to hydrogen ratio, atmospheric microturbulent velocity, mass- 
loss rate, and a measure of the acceleration of the transonic 
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outflow. Rotational velocities and terminal outflow velocities 
can be determined to considerable accuracy by means of ex- 
ternal methods such as rotational (de-) convolution methods 
(e.g.lHow arth et a ll 19971) and S EI-fitting of P-Cygni lines (e.g. 
Groenewegen & Lamers ; 1989), respectively. To get a good 
spectral fit it typically requires tens, sometimes hundreds of 
models per individual star. 

In the last few years the field of massive stars has seen 
the fortunate development that the number of O-type stars that 
have been studied spectroscopically has been doubled (e.g. 
Crowther etaD 12002 iHerrero et ail 120021 [Bianchi & Garcial 
2002| iBouret et alJl2003l iHillier et alJl2003| iGarcia & Bianchil 
2004; Martin s et al.l l2004 iMassev et all 120041 lEvans et alJ 
2004$. The available data set of massive O- and early B-type 
stars has recently again been doubled, mainly through the ad- 
vent of multi-object spectroscopy. Here we exp licitly men- 
tion t he VLT-FLAMES Survey of Massive Stars fevans et alJ 
120051) comprising over 100 hours of VLT time. In this survey 
multi-object spectroscopy using the Fibre Large Array Multi- 
Element Spectrograph (FLAMES) has been used to secure over 
550 spectra (of which in excess of 50 are spectral type O) in 
a total of seven clusters distributed over the Galaxy and the 
Magellanic Clouds. 

This brings within reach different types of studies that so 
far could only be attempted with a troublingly small sample 
of stars. These studies include establishing the mass loss be- 
haviour of Galactic stars across the upper Hertzsprung-Russell 
diagram, from the weak winds of the late O-type dwarfs (of 
order 10~ 8 Mgyr -1 ) to the very strong winds of early O-type 
supergiants (of order 10~ 5 M0yr~'); determination of the mass- 
loss versus metallicity dependence in the abundance range 
spanned by Small Magellanic Cloud to Galactic stars; placing 
constraints on the theory of massive star evolution by compar- 
ing spectroscopic mass determinations and abundance patterns 
with those predicted by stellar evolution computations, and the 
study of (projected) spatial gradients in the mass function of 
O- and B-type stars in young clusters, as well as such spatial 
gradients in the initial atmospheric composition of these stars. 

To best perform studies such as listed above not only re- 
quires a large set of young massive stars, it also calls for 
a robust, homogeneous and objective means to analyse such 
datasets using models that include state-of-the-art physics. This 
essentially requires an automated fitting method. Such an auto- 
mated method should not only be fast, it must also be suffi- 
ciently flexible to be able to treat early-type stars with widely 
different properties (e.g. mass-loss rates that differ by a factor 
of 10 3 ). Moreover, it should apply a well defined fitting cri- 
terium, like a x 2 criterium, allowing it to work in an automated 
and reproducible way. 

To cope with the dataset provided by the VLT-FLAMES 
Survey and to improve the objectivity of the analysis, we 
have investigated the possibility of automated fitting. Here we 
present a robust, fast, and accurate method to perform auto- 
mated fitting of the continuum normalized spectra of O- and 
early B-type stars with stellar winds using the fast pe rformance 
stellar atmosphere code fastwind fruls et alJl20"05l) combined 
with a genetic algorithm based fitting method. This first imple- 
mentation of an automated method should therefore be seen as 



an improvement over the standard "by eye" method, and not as 
a replacement of this method. The improvement lies in the fact 
that with the automated method large data sets (tens or more 
stars), spanning a wide parameter space, can be analysed in a 
repeatable and homogeneous way. It does not replace the "by 
eye" method as our automated fitting method still requires a 
by eye continuum normalization as well as a human controlled 
line selection. This latter should address the identification and 
exclusion of lines that are not modeled (i.e. blends), as well as 
introduce information on lacking physics and/or possible or po- 
tential problems in the model atmosphere code. Future imple- 
mentations of an automated fitting method may use the abso- 
lute spectrum, preferably over a broad wavelength range. This 
would eliminate the continuum rectification problem, however, 
it will require a modeling of the interstellar extinction. In this 
way one can work towards a true replacing of the "by eye" 
method by an automated approach. 

In Sect. |2 we describe the genetic algorithm method and 
implementation, and we provide a short resume of the applied 
unified, non-LTE, line-blanketed atmosphere code fastwind - 
which is the only code to date for which the method described 
here is actually achievable (in the context of analysing large 
data sets). To test the method we analyse a set of 12 early type 
spectra in Sect. [3] We start with a re-analysis of a set of seven 
stars in the open cl uster Cyg OB2 that have been studied by 
Herrero et al. (20 02). The advantage of focusing on this cluster 
is that it has been analysed with a previous version of fastwind, 
allowing for as meaningful a comparison as is possible, while 
still satisfying our preference to present a state-of-the-art anal- 
ysis. The analysis of Cyg OB2 has the added advantage that all 
stars studied are approximately equidistant. To test the perfor- 
mance of our method outside the parameter range offered by 
the Cyg OB2 sample we have included an additional five well- 
studied stars with either low density winds and/or very high 
rotational velocities. In Sect. 0] we describe our error analy- 
sis method for the multidimensional spectral fits obtained with 
the automated method. A systematic comparison of the ob- 
tained parameters with previously determined values is given 
in Sect. |5] Implications of the newly obtained parameters on 
the properties of massive stars are discussed in Sect. [6] In the 
last section we give our conclusions. 

2. Automated fitting using a genetic algorithm 

2.1. Spectral line fitting as an optimization problem 

Spectral line fitting of early-type stars is an optimization prob- 
lem in the sense that one tries to maximize the correspondence 
between a given observed spectrum and a synthetic spectrum 
produced by a stellar atmosphere model. Formally speaking 
one searches for the global optimum, i.e. best fit, in the param- 
eter space spanned by the free parameters of the stellar atmo- 
sphere model by minimizing the differences between the ob- 
served and synthesized line profiles. 

Until now the preferred method to achieve this minimiza- 
tion has been the so called fitting "by eye" method. In this 
method the best fit to the observed spectrum of a certain ob- 
ject is determined in an iterative manner. Starting with a first 
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guess for the model parameters a spectrum is synthesized. The 
quality of the fit to the observed spectrum is determined, as 
is obvious from the methods name, by an inspection by eye. 
Based on what the person performing the fit sees, for instance, 
whether the width of the line profiles are reproduced correctly, 
combined with his/her experience and knowledge of the model 
and the object, the model parameters are modified and a new 
spectrum is synthesized. This procedure is repeated until the 
quality of the fit determined by eye cannot be increased any- 
more by modifying the model parameters. 

It can be questioned whether a fit constructed with the fit- 
ting "by eye" method corresponds to the best fit possible, i.e. 
the global optimum. Reasons for this are, i) the restricted size 
of parameter space that can be investigated, both in terms of 
number of free parameters as well as absolute size of the pa- 
rameter domain that can investigated with high accuracy, ii) the 
limited number of free parameters that are changed simultane- 
ously, and Hi) biases introduced by judging the quality of a line 
fit by eye. The importance of the first point lies in the fact that 
in order to assure that the global optimum is found, a param- 
eter space that is as large as possible should be explored with 
the same accuracy for all parameters in the complete parame- 
ter space. If this is not the case the solution found will likely 
correspond to a local optimum. 

The argument above becomes stronger in view of the sec- 
ond point. Spectral fitting is a multidimensional problem in 
which the line profile shapes depend on all free parameters si- 
multaneously, though to a different extent. Consequently, the 
global optimum can only be found if all parameters are al- 
lowed to vary at the same time. The use of fit diagrams (e.g. 
iKudritzki & SimorJl978tlHerrero et all 19921) does not resolve 
this issue. These diagrams usually only take variations in T e g 
and \ogg into account, neglecting the effects of othe r pa- 
rameters, like microturbule nce (e.g. ISmith & Howarthl ll998: 
Villamariz & Herrerol HoOO ) and mass loss (e.g. [Fig. 5 of] 
Mokiem et all2004 . on the line profiles. 

The last point implies that, strictly speaking, fitting "by 
eye" cannot work in a reproducible way. There is no uniform 
well defined method to judge how well a synthetic line pro- 
file fits the data by eye. More importantly, it implies that there 
is no guarantee that the synthetic line profiles selected by the 
eye, correspond to the profiles which match the data the best. 
This predominantly increases the uncertainty in the derivation 
of those parameters that very sensitively react to the line profile 
shape, like for instance the surface gravity. 

The new fitting method presented here does not suffer from 
the drawbacks discussed above. It is an automated method ca- 
pable of global optimization in a multi-dimensional parameter 
space of arbitrary size ("Sect. 12.5V As it is automated, it does not 
require any human intervention in finding the best fit, avoiding 
potential biases introduced by "by eye" interpretations of line 
profiles. The method described here consists of two main com- 
ponents. The first component is the non-LTE stellar atmosphere 
code fastwind. Section l2~31 gives an overview of the capabilities 
of the code and the assumptions involved. The second compo- 
nent is the genetic algori t hm (G A) based optimizing routine 
pikaia from Charbonneau (1995), which is responsible for op- 
timizing the parameters of the fastwind models. For the tech- 



nical details of this routine and more information on GAs we 
refer to the cited paper and references therein. Here we will suf- 
fice with a short description of GAs and a description of the GA 
implementation with respect to optimization of spectral fits. 

2.2. The genetic algorithm implementation 

Genetic algorithms represent a class of heuristic optimization 
techniques, which are inspir ed by the notio n of evolution by 
means of natural selection (Darwin 1859). They provide a 
method of solving optimization problems by incorporating this 
biological notion in a numerical fashion. This is achieved by 
evolving the global solution over subsequent generations start- 
ing from a set of randomly guessed initial solutions, so called 
individuals. Selection pressure is imposed in between genera- 
tions based on the quality of the solutions, their so called fit- 
ness. A higher fitness implies a higher probability the solution 
will be selected for reproduction. Consequently, only a selected 
set of individuals will pass on their "genetic material" to sub- 
sequent new generations. 

To create the new generations discussed above GAs require 
a reproduction mechanism. In its most basic form this mecha- 
nism consists of two genetic operators. These are the crossover 
operator, simulating sexual reproduction, and the mutation op- 
erator, simulating copying errors and random effects affecting 
a gene in isolation. An important benefit of these two opera- 
tors is the fact that they also introduce new genetic material 
into the population. This allows the GA to explore new re- 
gions of parameters space, which is important in view of the 
existence of local extremes. When the optimization runs into 
a local optimum, these two operators, where usually muta- 
tion has the strongest effect, allow for the construction of in- 
dividuals outside of this optimum, thereby allowing it to find 
a path out of the local optimum. This capability to escape lo- 
cal extremes, consequently, classifies GAs as global optimiz- 
ers and is one of the reasons they have been applied to many 
problems in and outside astrop hysics (e.g. Metcal fe et al . 2000; 
IGibson & Charbonneaull998l) . 

Using an example we can further illustrate the GA opti- 
mization technique. Lets assume that the optimization problem 
is the minimization of some function /. This function has n 
variables, serving as the genetic building blocks, spanning a 
n dimensional parameter space. The first step in solving this 
problem is to create an initial population of individuals, which 
are sets of n parameters, randomly distributed in parameter 
space. For each of these individuals the quality of their solution 
is determined by simply calculating / for the specific param- 
eter values. Now selection pressure is imposed and the fittest 
individuals, i.e. those that correspond to the lowest values of /, 
are selected to construct a new generation. As the selected in- 
dividuals represent the fittest individuals from the population, 
every new generation will consist of fitter individuals, leading 
to a minimization of /, thereby, solving the optimization prob- 
lem. 

With the previous example in mind we can explain our im- 
plementation of the GA for solving the optimization problem 
of spectral line fitting, with the following scheme. We start out 
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with a first generation of a population of fastwind models ran- 
domly distributed in the free parameter space (see Sect. 12. 5> . 
For each of these models it is determined how well an observed 
spectrum is fitted by calculating the reduced chi squared, x] eA ■, 
for each of the fitted lines i. The fitness F, of a model is then 
defined as the inverted sum of the^ ed ,'s, i.e. 



with the number of processors used. Test calculations showed 
that for configurations in which the population size is an integer 
multiple of the number of processors the sequential overhead is 
negligible. Consequently, the runtime scales directly with the 
inverse of the number of processors. Thus enabling the auto- 
mated fitting of spectra. 
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(1) 2.4. The non-LTE model atmosphere code fastwind 



where N corresponds to the number of lines evaluated. The 
fittest models are selected and a new generation of models is 
constructed based on their parameters. From this generation the 
fitnesses of the models are determined and again from the fittest 
individuals a new generation is constructed. This is repeated 
until F is maximized, i.e. a good fit is obtained. 

In terms of quantifying the fit quality Eq. does not rep- 
resent a unique choice. Other expressions for the fitness cri- 
terium, for instance, the sum of the inverted x*. .'s of the indi- 
vidual lines, or the inverted of all the spectral points eval- 
uated, also produce the required functionality of an increased 
fitness with an increased fit quality. We have chosen this par- 
ticular form based on two of its properties. Firstly, the eval- 
uation of the fit quality of the lines enter into the expression 
individually, ensuring that, regardless of the number of points 
in a certain line, all lines are weighted equally. This allows as 
well for weighting factors for individual lines, which express 
the quality with which the stellar atmosphere synthesizes these 
lines (cf. Sect. 13. 2\ . Secondly, using the inverted sum of the 
Xte&' s msteac l °f me sum °f me mverte d^ ed ,'s avoids having 
a single line, which is fitted particularly well, to dominate the 
solution. Instead the former form demands a good fit of all lines 
simultaneously. 

2.3. Parallelization of the genetic algorithm 

The ability of global optimization of GAs comes at a price. 
Finding the global minimum requires the calculation of many 
generations. In Sect. l2.6l we will show that for the spectra stud- 
ied in this paper, the evaluation of more than a hundred genera- 
tions is needed to assure that the global optimum is found. For 
a typical population size of -70 individuals, this comes down 
to the calculation of ~7000 fastwind models. With a modern 
3 GHz processor a single fastwind model (aiming at the analy- 
sis of hydrogen and helium lines) can be calculated within five 
to ten minutes. Consequently, automated fitting on a sequential 
computer would be unworkable. 

To overcome this problem, parallelization of the pikaia rou- 
tin e is necessary. This p arallelization is inspired by the work 
of iMetcalfe & Charbonneaul J2003h . Consequently, our paral- 
lel version is very similar to the version of these authors. The 
main difference between the two versions, is an extra paral- 
lelization of the so called elitism opt ion in the reproduction 
schemes (see IMetcalfe & C harbonn eauj). This was treated in a 
sequential manner in the IMetcalfe & Charbonneaul implemen- 
tation and has now been parallelized as well. 

Due to the strong inherent parallelism of GAs, the paral- 
lel version of our automated fitting method scales very well 



For modeling the optical spectra of our stars we use the latest 
version of the non-LTE, line-blanketed atmosphere code fast- 
wind for e arly-type stars wit h winds. For a detailed description 
we refer to lPulsetalJ £2005 ). Here we give a short overview of 
the assumptions made in this method. The code has been devel- 
oped with the emphasis on a fast performance (hence its name), 
which makes it currently the best suited (and realistically only) 
model for use in this kind of automated fitting methods. 

fastwind adopts the concept of "unified model atmo- 
spheres", i.e. including both a pseudo-hydrostatic photosphere 
and a transonic stellar wind, assuring a smooth transition be- 
tween the two. The photospheric density structure follows from 
a self-consistent solution of the equation of hydrostatic equilib- 
rium and accounts for the actual temperature stratification and 
radiation pressure. The temperature calculation utilizes a flux- 
correction method in the lower atmosphere and the thermal bal- 
ance of electrons in the outer atmosphere (with a lower cut-off 
at r m i n = 0.5r e ff). In the photosphere the velocity structure, 
v(r), corresponds to quasi-hydrostatic equilibrium; outside of 
this regime, in the region of the sonic velocity and in the super- 
sonic wind regime it is prescribed by a standard /3-type velocity 
law, i.e. 



(2) 



where Voo is the terminal velocity of the wind. The parameter 
r is used to assure a smooth connection, and f3 is a measure of 
the flow acceleration. 

The code distinguishes between explicit elements (in our 
case hydrogen and helium) and background elements (most im- 
portantly: C, N, O, Ne, Mg, Si, S, Ar, Fe, Ni). The explicit el- 
ements are used as diagnostic tools and are treated with high 
precision, i.e. by detailed atomic models and by means of co- 
moving -frame transport for the line transitions. The Hi and 
He ii model atoms consist of 20 levels each; the He i model in- 
cludes levels up to and including n = 10, where levels with 
n > 8 have been packed. The background ions are included to 
allow for the effects of line-blocking (treated in an approximate 
way by using suitable means for the corresponding line opaci- 
ties) and line-blanketing. Occupation numbers and opacities of 
both the explicit and the most abundant background ions are 
constrained by assuming statistical equilibrium. The only dif- 
ference between the treatment of these types of ions is that for 
the background ions the Sobolev approximation is used in de- 
scribing the line transfer (accounting for the actual illumination 
radiation field). 

Abundances of the ba ckgroun d elements a re take n from the 
solar values provided bv lGrevesse & Sauvall (11998. and refer- 
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ences therein). The He/H ratio is not fixed and can be scaled 
independently from the background element abundances. 

A comparison between the optical H and He lines as syn- 
thesized by fastwind and those predicted by the independent 
comparison code cmfgen jHillier & MilleJfl998h show excel- 
lent agreement, save for the He i singlet lines in the temperature 
range between 36 000 and 41 000 K for dwarfs and between 
3 1 000 and 35 000 K for supergiants, where cmfgen predicts 
weaker lines. We give account of this discrepancy, and there- 
fore of an increased uncertainty in the reproduction of these 
lines, by introducing weighting factors, which for the He i sin- 
glets of stars in these ranges are lower (cf. Sect. 13. 2> . 

2.5. Fit parameters 

The main parameters which will be determined from a spectral 
fit using fastwind are the effective temperature T e ff, the surface 
gravity g, the microturbulent velocity v tul -b, the helium over hy- 
drogen number density Yn e , the mass loss rate M and the expo- 
nent of the beta-type velocity law /?. These parameters span the 
free parameter space of our fitting method. The stellar radius, 
R+, is not a free parameter as its value is constrained by the 
absolute visual magn itude My. To cal culate we adopt the 
procedure outlined in iKudritzkil Jl9 80). i.e. 

5 log R/Rq = 29.57 - (M v - V) , (3) 

where V is the visual flux of the theoretical model given by 

-2.5 log F A S A dA. (4) 
Jo 

In the above equa t ion S * is the V-filter function of 
Matt hews & Sandag el (fl963) and Fx is the theoretical stellar 
flux. Note that as R* is an input parameter, Fx is not known 
before the fastwind model is calculated. Therefore, during the 
automated fitting w e approximate F^ by a black body radiating 
at T = 0.9r eff (cf. iMarkova et alll2004l) . After the fit is com- 
pleted we use the theoretical flux from the best fit model to cal- 
culate the non approximated stellar radius. Based on this radius 
we rescale th e mass loss rate using the invariant wi nd-strength 
parameter Q dPuls et alJl996Hde Koter et alll997t) 



The largest difference between the approximated and final stel- 
lar radius for the objects studied here, is ~2 percent. The corre- 
sponding rescaling in M is approximately three percent. 

The projected rotation velocity, v r sin i, and terminal veloc- 
ity of the wind are not treated as free parameters. The value of 
v r sin i is determined from the broadening of weak metal lines 
and the width of the He i lines. For v ro we adopt values ob- 
tained from the study of ultraviolet (UV) resonance lines, or, if 
not available, values from calibrations are used. 

Our fitting method only requires the size of the free param- 
eter domain to be specified. For the objects studied in this paper 
we keep the boundaries between which the parameters are al- 
lowed to vary, fixed for vtm-b, Y^ e and /3. The adopted ranges, 
respectively, are [0, 20] kms" 1 , [0.05, 0.30] and [0.5, 1.5]. The 



boundaries for T e ff are set based on the spectral type and lumi- 
nosity class of the studied object. Usually the size of this range 
is set to approximately 5000 K. The logg range is delimited 
so that the implied stellar mass lies between reasonable bound- 
aries. For instance for the B 1 1 star Cyg OB2 #2 the adopted T e g 
range together with its absolute visual magnitude imply a pos- 
sible range in/?* of [1 1.5:12. Q]Rq. For the automated fit we set 
the minimum and maximum logg to 3.1 and 3.8, respectively, 
which sets the corresponding mass range that will be investi- 
gated to [5.0:25.2] Mq. For the mass loss rate we adopt a con- 
servative range of at least one order of magnitude. As example 
for the analysis of Cyg OB2 #2 we adopted lower and upper 
boundaries of 4 x 10~ 8 and 2 x 10~ 6 Mgyr -1 , respectively. 

2.6. Formal tests of convergence 

Before we apply our automated fitting method to real spec- 
tra, we first test whether the method is capable of global opti- 
mization. For this we perform convergence tests using synthetic 
data. The main goal of these tests is to determine how well and 
how fast the input parameters, used to create the synthetic data, 
can be recovered with the method. The speed with which the 
input parameters are recovered, i.e. the number of generations 
needed to find the global optimum, can then be used to deter- 
mine how many generations are needed to obtain the best fit for 
a real spectrum. In other words, when the fit has converged to 
the global optimum. 

Three synthetic datasets, denoted by A, B and C, were cre- 
ated with the following procedure. First, line profiles of Balmer 
hydrogen lines and helium lines in the optical blue and Ha 
in the red calculated by fastwind were convolved with a ro- 
tational broadening profile. Tabled lists the parameters of the 
three sets of models as well as the projected rotational veloc- 
ity used. A second convolution with a Gaussian instrumental 
profile was applied to obtain a spectral resolution of 0.8 A and 
1.3 A for, respectively, the Ha line and all other lines. These 
values correspond to the minimum resolution of the spectra fit- 
ted in Sect. [3] Finally, Gaussian distributed noise, correspond- 
ing to a signal to noise value of 100, was added to the profiles. 
Dataset A represents an 03 I star with a very dense stellar wind 
(M = 10~ 5 M Q yr _1 ), while set B is that of an 05.5 I with a 
more typical O-star mass loss. The last set C is characteristic 
for a B0 V star with a very tenuous wind of only 10~ 8 MQyr -1 . 

From the synthetic datasets we fitted nine lines, three hy- 
drogen, three neutral helium and three singly ionized helium 
lines, corresponding to the minimum set of lines fitted for a 
single object in Sect. [3] The fits were obtained by evolving a 
population of 72 fastwind models over a course of 200 genera- 
tions. In this test and throughout the remainder of the paper we 
use pikaia with a dynamically adjustable mutation rate, with the 
minimum and maximum mu tation rate set to the default values 
(see Charbonneau & Knapp 1995). Selection pressure, i.e. the 
weighting of the probability an individual will be selected for 
reproduction based on its fitness, was also set to the default 
value. 

Table^lists the parameter ranges in which the method was 
allowed to search, i.e. the minimum and maximum values al- 



6 



M. R. Mokiem et al.: Spectral analysis using a genetic algorithm based fitting method 



Table 1. Input parameters of the formal test models ("In" column) and parameters obtained with the automated fitting method by 
fitting synthetic data created from these models ("Out" column). Results were obtained by evolving a population of 72 fastwind 
models over 200 generations. 





Set A 


Search 




SetB 


Search 




SetC 


Search 






In 


range 


Out 


In 


range 


Out 


In 


range 


Out 


Spectral type 


031 






05.51 






B0 V 






r cff [kK] 


45.0 


[42, 47] 


45.0 


37.5 


[35, 40] 


37.6 


30.0 


[28, 34] 


29.9 


logg [cm s~ 2 ] 


3.80 


[3.5, 4.0] 


3.84 


3.60 


[3.3, 3.9] 


3.57 


4.00 


[3.7, 4.3] 


3.95 




17.0 






20.0 






8.0 






logi* [Lq] 


6.03 






5.85 






4.67 






Vturb [kms" 1 ] 


5.0 


[0, 20] 


5.9 


10.0 


[0, 20] 


9.7 


15.0 


[0, 20] 


14.8 




0.15 


[0.05, 0.30] 


0.15 


0.10 


[0.05, 0.30] 


0.10 


0.10 


[0.05, 0.30] 


0.10 


M [lO^Moyr 1 ] 


10.0 


[1.0, 20.0] 


9.3 


5.0 


[1.0, 10.0] 


5.3 


0.01 


[0.001, 0.2] 


0.008 


P 


1.20 


[0.5, 1.5] 


1.18 


1.00 


[0.5, 1.5] 


0.99 


0.80 


[0.5, 1.5] 


0.93 


v„ [kms -1 ] 


2500 






2200 






2000 






v, sin i [km s~'] 


150 






120 






90 







lowed for the parameters of the fastwind models. As Voo and 
v r sin i are not free parameters these were set equal to the input 
values. 

In all the three test cases the automated method was able to 
recover the global optimum. TableJ^lists the parameters of the 
best fit models obtained by the method in the "Out" columns. 
Compared to the parameters used to create the synthetic data, 
there is very good agreement. Moderate differences (of a 15- 
20% level) are found for v tur b recovered from dataset A and 
for the wind parameters f3 and M recovered from dataset C. 
This was to be expected. In the case of the wind parameters 
the precision with which information about these parameters 
can be recovered from t he line profiles decreases with decreas- 
ing wind density (e.g. Pu is et al . 1996). Still, the precision with 
which the wind parameters are recovered for the weak wind 
data set C, is remarkable. 

A similar reasoning applies for the microturbulent veloc- 
ity recovered from data set A. For low values of the micro- 
turbulence, i.e. Vturb < v t h, thermal broadening will dominate 
over broadening due to microturbulence. This decreases the 
precision with which this parameter can be recovered from the 
line profiles. Realizing that in case of this dataset for helium 
v t h ~ 14 km s _1 , again, the precision with which Vturb is recov- 
ered, is impressive. 

To illustrate how quickly and how well the input parame- 
ters are recovered Fig.[^shows the evolution of the fit parame- 
ters during the fit of synthetic dataset B. Also shown, as a grey 
dashed line, is the fitness of the best fitting model found, dur- 
ing the run. This fitness is normalized with respect to the fitness 
of the model used to create the synthetic data (the data being 
the combination of this model and noise). Note that the final 
maximum normalized fitness found by the method exceeds 1 .0, 
which is due to the added noise allowing a further fine tuning 
of the parameters by the GA based optimization. As can be 



seen in this figure the method modifies multiple parameters si- 
multaneously to produce a better fit. This allows for an efficient 
exploration of parameter space and, more importantly, it allows 
for the method to actually find the global optimum. 

In the case of dataset B finding the global optimum required 
only a few tens of generations (~30). For the other two datasets 
all save one parameter were well established within this num- 
ber of generations. To establish the very low value of v tur b in 
dataset A and the very low M in dataset C required -100 gen- 
erations. We will adopt 150 generations to fit the spectra in 
Sect. One reason, obviously, is that to safeguard that the 
global optimum is found. A second reason, however, is that it 
assures that the errors on the model parameters that we deter- 
mine are meaningful (i.e. it assures that the error on the error is 
modest). 

We consider doing such a formal test as performed above 
as part of the analysis of a set of observed spectra, as the exact 
number of generations required is, in principle, a function of 
e.g. the signal-to-noise ratio and the spectral resolution. Also, 
special circumstances may play a role, such as potential nebular 
contamination (in which case the impact of removing the line 
cores from the fit procedure needs to be assessed). 

3. Spectral analysis of early-type stars 

In this section we apply our fitting method to sev en stars in the 
open c luster Cyg OB2, previously analysed by iHerrero et alJ 
(2002) and five "standard" early-type stars, 10 Lac, t Sco, 
£ Oph, HD 15629 and HD 217086, previously analysed by var- 
ious authors. 

3.1. Description of the data 

Table|2lists the basic properties of the data used for the analy- 
sis. All spectra studied have a S/N of at least 100. The spectral 
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generation generation 



Fig. 1. Evolution of the best fitting model parameters for formal test B. From the 200 generation run only the first 75 generations 
are shown. For this specific data set the location of the global optimum is found within 50 generations. This is indicated by the 
highest fitness found during the run, which is shown as a grey dashed line and is scaled to the right vertical axis. The fitness is 
normalized with respect to the fitness of the model used to create the synthetic data (the data being this model plus noise). 



resolution of the data in the blue (regions between ~4000 and 
~5000 A) and the red (region around Ha) is given in Tab.|2] 

The optical spectra o f the stars in Cyg OB2 we re obtained 
by IHerrero et ail {l999) and iHerrero et all ll2000l) . Absolute 
visual magnit udes of the Cyg OB2 objects were adopted 
from Massev & Thompson (1991), and correspond to a dis- 
tance modulus of 11 . 2 m . Note that for object #8A Tab. 7 in 
iMassev & Tho mpson contains an incorrect Vn value of 4.08 m . 
This should have been 4.26 m conform the absorption given in 
this table and the visual magnitude i n thei r Tab. 2. For v, sin / 
values determined by IHerrero et al J {.2002) are used, with the 
exception of objects #8A and #10. For these we found that 
the He i and metal lines are somewhat better reproduced if 
we adopt v, -sin / that are higher by ~35% and ~10%, respec- 
tively. Terminal flow velocities of the wind have been obtained 
from UV spectra o btained with Hubble Space Telescope (cf. 
IHerrero et alfcOOlh. Data of H P 1 5629, HP 2170 86 and ( Oph 
are from Herr ero et alJ 1 1992) and Herrero ( 1993). For My, 
and v r sin / values given bv |Repolust et al. (2004) are adopted. 



The distances to these objects are based on spectroscopic paral- 
laxes, except for C O ph which has a reliable Hipparcos distance 
JSchroderetalJ2004l). 

Th e spectrum of 10 Lac was obtained by IHerrero et al.l 
J2002I) . The absolute visual magnitude of this star is from 



Herre ro et al 1 dl992l) . For Voo we adopted the minimum value 
which is approximately equal to the escape velocity at the stel- 
lar surface of this object. For the projected rotational v elocity 
we ad opt 35 km s -1 . The blue spectrum of t Sco is f romlKiliai] 
dl992l) . The red region around Ha was observed bv lZaal et alJ 
( 1999). For r Sco we also adopt the Hipparcos distance. This 
distance results in an absolute visual magnitude which is rather 
large for the sp ectral type of this object, but is in bet ween the 
My adopted by Kilian ( 1993) and iHumphrevsl 11978). For the 
projected rotational velocity a value of 5 km s _1 was adopted. 

3.2. Lines selected for fitting and weighting scheme 

For the analysis fastwind will fit the hydrogen and helium spec- 
trum of the investigated objects. Pepending on the wavelength 
range of the available data, these lines comprise for hydrogen 
the Balmer lines Ha, H/3, Hy and HS; for He i the singlet lines at 
4387 and 4922 A, the He i triplet lines at 4026, which is blended 
with He ii, 4471 and 4713 A; and finally for Hen the lines at 
4200, 4541 and 4686 A. 

For an efficient and reliable use of the automated method 
we have to incorporate into it the expertise that we have devel- 
oped in the analysis of OB stars. The method has to take into 
account that some lines may be blended or that they cannot 
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Table 2. Basic param eters of the early type stars s tudied here. Spectral types are taken from iMassev & Thompsonl Jl99ll) . 
Walbornl (ll972lll973l) and IConti & Alschulerl dl97 lh . Blue and red resolution, respectively, correspond to the region between 
~4000 and ~5000 A and the region around Ha. 



Star 


Spectral 


M v 


Blue 


Red 


v r sin i 


Voo 




Type 




resolution [A] 


resolution [A] 


[kms- 1 ] 


[kms- 1 ] 


Cyg OB2 #7 


03 If 


-5.91 


0.6 


0.8 


105 


3080 


CygOB2#ll 


05 If + 


-6.51 


1.3 


0.8 


120 


2300 


Cyg OB2 #8C 


05 If 


-5.61 


1.3 


0.8 


145 


2650 


Cyg OB2 #8A 


05.5 1(f) 


-6.91 


0.6 


0.8 


130 


2650 


Cyg OB2 #4 


07 ni((f)) 


-5.44 


1.3 


0.8 


125 


2550 


Cyg OB2 #10 


09.51 


-6.86 


0.6 


0.8 


95 


1650 


Cyg OB2 #2 


Bl I 


-4.64 


0.6 


0.8 


50 


1250 


HD 15629 


05 V((f)) 


-5.50 


0.6 


0.8 


90 


3200 


HD 217086 


07 Vn 


-4.50 


0.6 


0.8 


350 


2550 


10 Lac 


09 V 


-4.40 


0.6 


0.6 


35 


1140 


^Oph 


09 V 


-4.35 


0.6 


0.8 


400 


1550 


T SCO 


B0.2 V 


-3.10 


0.2 


0.2 


5 


2000 



Table 3. Line weighting scheme adopted for different spectral 
types and luminosity classes for the objects fitted in this paper. 
Late, mid and early spectral type correspond to, respectively, 
[02-05.5], [06-07.5] and [08-B1]. The weights are imple- 
mented in the fitness definition according to Eq. (Jfji and have 
values of 1.0, 0.5 and 0.25 in case of h, m and 1, respectively. 



have assigned the spectral lines different weights depending 
on their behaviour in each stellar group. This behaviour rep- 
resents the expertise from years of "by eye" data analysis that 
is being translated to the method. Three different weights are 
assigned to each line: high, to lines very reliable for the anal- 
ysis; medium, and low. The implementation of these weights 
into the fitness definition is given by 



Dwarfs 
Late Mid Early 



Super Giants 
Late Mid Early 



F = 



H B aimer 


h 


h 


h 


h 


h 


h 


He i singlets 


h 


1 


1 


h 


1 


1 


He i 4026 


h 


h 


h 


h 


h 


h 


He i 4471 


h 


h 


h 


1 


m 


h 


Hei4713 


h 


h 


h 


h 


h 


h 


Hen 4686 


h 


m 


m 


m 


m 


m 


Hen 4541 


h 


h 


h 


h 


h 


h 


Hen 4200 


m 


m 


m 


m 


m 


m 



be completely reproduced by the model atmosphere code for 
whatever r eason For ex ample, the so-called "generalized dilu- 
tion effecf' IVoels et alJ lll989l) . present in the He i A441 1 line in 
late type supergiants, that is still lacking an explanation. 

To that end we have divided the stars in two classes 
("dwarfs" and "supergiants", following their luminosity class 
classification 1 ), and three groups in each class (following spec- 
tral types). We have then a total of six stellar groups, and 



1 For the one giant in our sample, Cyg OB2 #4, we have adopted 
the line weighting scheme for dwarfs. 



(6) 



where the parameter w, corresponds to the weight of a specific 
line. 

Table|2gives the weights assigned to each line in each stel- 
lar group. We will only briefly comment on the low or medium 
weights. He i singlets are assigned a low weight for mid-type 
stars because of the singlet differentia l beha viour found be- 
tween fastwind and cmfgen dPuls et al.ll2005t) . while they are 
very weak for early-type stars. In these two cases therefore we 
prefer to rely on the triplet He 1/14471 line. To this line, how- 
ever, a low weight is assigned at late-type Supergiants because 
of the above mentioned dilution effect. 

He ii /14686 is only assigned a medium weight (except for 
late type dwarfs), as this line is not always completely consis- 
tent with the mass-loss rates derived from Halpha. He n ,44200 
is sometimes blended with N m /14200, and sometimes it is not 
completely consistent with the rest of the He n lines. He i and 
He ii lines at 4026 A do overlap, but for both lines we find a 
consistent behaviour. 

The highest weight is therefore given to the Balmer lines 
plus the He n /14541 and the He i/He n 4026 lines, which define 
the He ionization balance with He 1/14471 or the singlet Hei 
lines. Note however that, as discussed above, all lines fit simul- 
taneously in a satisfactory way for our best fitting models. 
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Table 4. Results obtained for the investigated early type stars using GA optimized spectral fits. The spectra were fitted by 
evolving a population of 72 fastwind models over a course of 150 generations. Spectros copic masses M 5 are calculated with the 
gravities corrected for centrifugal acceleration logg c . Evolutionary masses M ev are from lSchaller et all i 19921) . The error bars on 
the derived parameters are given in Tab.|5]and are discussed in Sect.@] 



Star 


T e s 


logg 


log ft 


R* 


log L* 




Vturb 


M 




P 


M, 


M ev 




[kK] 


[cm s ] 


[cm s" 2 ] 


[Rq] 


[Lq] 




[kms 1 ] 


[M©yr 


J 




[Mq] 


[Mq] 


Cyg OB2 #7 


45.8 


3.93 


3.94 


14.4 


5.91 


0.21 


19.9 


9.98-10- 




0.77 


65.1 


67.8 


Cyg OB2#ll 


36.5 


3.62 


3.63 


22.1 


5.89 


0.10 


19.8 


7.36-10- 


-6 


1.03 


75.9 


55.6 


Cyg OB2 #8C 


41.8 


3.73 


3.74 


13.3 


5.69 


0.13 


0.5 


3.37-10- 


-6 


0.85 


36.0 


49.2 


Cyg OB2 #8A 


38.2 


3.56 


3.57 


25.6 


6.10 


0.14 


18.3 


1.04-10" 


-5 


0.74 


89.0 


74.4 


Cyg OB2 #4 


34.9 


3.50 


3.52 


13.7 


5.40 


0.10 


18.9 


8.39-10" 


-7 


1.16 


22.4 


32.5 


Cyg OB2 #10 


29.7 


3.23 


3.24 


29.9 


5.79 


0.08 


17.0 


2.63-10" 


-6 


1.05 


56.0 


45.9 


Cyg OB2 #2 


28.7 


3.56 


3.57 


11.3 


4.88 


0.08 


16.5 


1.63-10 


-7 


0.80" 


17.0 


18.7 


HD 15629 


42.0 


3.81 


3.82 


12.6 


5.64 


0.10 


8.6 


9.28-10" 


-7 


1.18 


37.8 


47.4 


HD 217086 


38.1 


3.91 


4.01 


8.30 


5.11 


0.09 


17.1 


2.09-10" 


-7 


1.27 


25.7 


28.5 


10 Lac 


36.0 


4.03 


4.03 


8.27 


5.01 


0.09 


15.5 


6.06-10" 


-8 


0.80" 


26.9 


24.9 


<TOph 


32.1 


3.62 


3.83 


8.9 


4.88 


0.11 


19.7 


1.43-10 


-7 


0.80" 


19.5 


20.3 


T SCO 


31.9 


4.15 


4.15 


5.2 


4.39 


0.12 


10.8 


6.14-10" 


-8 


0.80" 


13.7 


16.0 



" assumed fixed value 



3.3. Fits and comments on the individual analysis 

In the following we will present the fits that were obtained by 
the automated method for our sample of 12 early type stars, 
and comment on the individual analysis of the objects. Listed 
in Tab. 0] are the values determined for the six free parameters 
investigated and quantities derived from these. 



3.3.1 . Analysis of the Cyg OB2 stars 

Th e Cyg OB2 objects s tudied he re wer e previously analysed 
bv lHerreroetall ( l2002l hereafter IhPnT) . We opted to reanal- 
yse these stars (to test our method) as these stars have equal 
distances and have been analysed in a homogeneous way us- 
ing (an earlier version of) the same model atmosphere code. In 
Sect. |3 we will sy stematically compare our results with those 
obtained bv lHPNl Here, we will incidentally discuss the agree- 
ment if this turns out to be relatively poor or if the absolute 
value of a parameter seems unexpected, and we wanted to test 
possible causes for the discrepancy. 

Cyg OB2 #7 The best fit obtained with our automated fitting 
method for Cyg OB2 #7 is shown in Fig. [2] For all hydrogen 
lines fitted, including HS not shown here, and all He n lines the 
fits are of very good quality. Note that given the noise level the 
fits of the He i lines are also acceptable. 

Interesting to mention is the manner in which the He i and 
Hen blend at 4026 A is fitted. At first sight, i.e. "by eye", it 
seems that the fit is of poor quality, as the line wings of the 
synthetic profile runs through "features" which might be at- 
tributed to blends of weak photospheric metal lines. However, 
the broadest of these features have a half maximum width of 



~70 kms , which is much smaller than the projected rota- 
tional velocity of 105 kms -1 . Consequently, these features are 
dominated by pure noise. 

Compared to the investigation of HPN we have partial 
agreement between the derived parameters. The mass loss rate, 
Jeff and to a lesser degree /3 agree very well. For \ogg and 
the helium abundance we find, however, large differences. The 
log g value obtained here is ~0.2 dex larger, which results in 
a spectroscopic mass of 65.1 Mq. A value which is in good 
agreement with the evolutionary mass of 67.8 Mq. 

The helium abundance needed to fit this object is 0.21 , 
which is considerably lower than the value obtained bv IHPNI 
who found an abundance ratio of 0.31. This large value still 
corresponds to a strong helium surface enrichment. An inter- 
esting question we need to address, is whether this is a real 
enrichment and not an artifact that is attributable to a degen- 
eracy effect of T e s and Yn e . The latter can be the case, as no 
He i lines are present in the optical spectrum of Cyg OB2 #7. 
This issue can be resolved with our fitting method by refitting 
the spectrum with a helium abundance fixed at a lower value 
than previously obtained. If T^g and Yu e are truly degenerate 
this would again yield a good fit, however for a different T e ff. 

Shown as dotted lines in Fig. |2 are the results of refitting 
Cyg OB2 #7 with a helium abundance fixed at the solar value. 
For this lower Yu e a 7" e ff that is lower by ~2. 1 kK was obtained. 
This was to be expected as for this temperature regime He m 
is the dominant ionization stage. When consequently Fne is re- 
duced a reduction of the temperature is required to fit the He n 
lines. The reduction of T e g obtained is the maximum for which 
still a good fit of the hydrogen lines is possible and the Hei 
lines do not become too strong. More importantly, in Fig. [2] it 
is shown that even with this large reduction of T e g the Hen 
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Hel 4471 Hel 4922 

1.02 [ \ 1.02 




6520 6560 6600 
Hell 4200 




0.98 

4020 4025 4030 4384 4388 4392 

Hell 4541 Hell 4686 

1.15 



1.10 





4915 4920 4925 



4190 4200 4210 
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Fig. 2. Comparison of the observed line profiles of Cyg OB2 #7 with the best fit obtained by the automated fitting method 
(dashed lines). Note that the He n line at 6527.1 A is not included in the fit and, therefore, disregarded by the automated method. 
Horizontal axis gives the wavelength in A. Vertical axises give the continuum normalized flux and are scaled differently for each 
line. In this figure the dotted lines correspond to a fit obtained for a helium abundance fixed at 0. 1 . See text for further comments. 
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Fig. 3. Same as Fig.|2 however for Cyg OB2 #1 1 



lines cannot be fitted. This implies that T e g and Fee are not 
degenerate and the obtained helium enrichment is real. 

Cyg OB2 #1 1 Figure|3]shows the fit to Cyg OB2 #1 1 . In gen- 
eral all lines are reproduced correctly. There is a slight un- 
der prediction of the cores of Hy and Hen A4541 , a problem 
that was also pointed out by iHerrero et alJ i 19921) and IhPNL 
Possibly this is due to too much filling in of the predicted 
profiles by wind emission. Part of the Hen /14541 discrep- 
ancy might be related to problems in the theoretical broadening 
functions (see lReoolust et al.l2005l) . 

The parameters obtained for this object, with exception of 
M, are in agreement with the parameters derived bv lHPNl With 
our automated method a mass loss rate lower by ~0. 1 dex was 
obtained. Note that due to this lower value the behaviour of this 



object in terms of its modified wind momentum (cf. Sect. 15. 4> 
is in better accord with that of the bulk of the stars investigated 
in this paper. 

Cyg OB2 #8C The best fit for Cyg OB2 #8C is shown in 
Fig. |3 Again, with exception of the mass loss rate, the param- 
eters we obtain for this object are in good agreement with the 
findings of lHPNl We do find a small helium abundance en- 
hancement, whereas HPN found a solar value. 

To fit the P Cygni type profile of He n at 4686 A, the au- 
tomated method used a M which, compared to these authors, 
was higher by approximately 0.15 dex. This higher value for 
the mass loss rate results in a Ha profile which, at first sight, 
looks to be filled in too much by wind emission. To assess 
whether this could correspond to a significant overestimation of 
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Fig. 4. Same as Fig. [2] however for Cyg OB2 #8C. Shown with a dotted line for Ha and Hen A4686 are the line profiles of a 
model with a 0.05 dex lower M, which "by eye" fits the core of Ha. See text for further comments. 



H5 Hy Hp Ha Hel 4387 




4090 4100 4110 4330 4340 4350 4850 4860 4870 6530 6560 6590 4384 4388 4392 

Hel 4471 Hel 4922 Hell 4200 Hell 4541 Hell 4686 
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Fig. 5. Same as Fig.|2l however for Cyg OB2 #8A. The dotted lines correspond to a model with a M higher by 0.04 dex. This 
mass loss rate was obtained by fitting the best fit model, found by the automated method, "by eye" to the Ha core. Even though 
the fit obtained with the higher M results in a fit of Ha which is more pleasing to the eye in the line core, this higher mass loss 
rate does not describe this object the best. This can be seen best from the reduced fit quality of the other hydrogen Balmer lines 
and the severe mismatch of He i ,4447 1 . See text for further comments. 



the mass loss rate, we lowered M in the best fit model by hand 
until the core of Ha was fitted. In Fig.@]the resulting line pro- 
files are shown as a dotted line for Ha and He n /14686, which 
for this fit are the lines which visibly reacted to the change in 
mass loss rate. To obtain this fit "by eye" of the Ha core, a 
reduction of M with merely 0.05 dex was required, showing 
that the mass loss rate was not overestimated by the automated 
method. Note that for this lower mass loss rate the fit of the 
He ii /14686 becomes significantly poorer. 

Cyg OB2 #8A be Becker etaD <2004 report this to be a 06 1 
and 05.5 III binary system, therefore the derived parameters, 
in particular the spectroscopically determined mass, should be 



taken with care. However, as this paper also aims to test auto- 
mated fitting we did pursue the comparison of this object with 
HPN, who also treated the system assuming it to be a single 
star. 

We obtained a good fit for all lines except for the problem- 
atic Hen /14686 line. The best fit is shown in Fig. [5] Again the 
Ha core is not fitted perfectly. To determine how significant this 
small discrepancy is, we fitted the Ha core in a similar manner 
as for Cyg OB2 #8C. To obtain a good fit "by eye" we find that 
M has to be increased by 0.04 dex, indicating the extreme sen- 
sitivity of Ha to M in this regime. The profiles corresponding 
to the increased mass loss rate model are shown in Fig. [5] as 
dotted lines. It is clear that not only the "classical" wind lines 
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Fig. 7. Same as Fig.|2] however for Cyg OB2 #10. The emission feature in the core of Ha was not included in the fit. A subsequent 
test which did include this feature in the fit yielded the same parameters except for a small increase of M with 0.04 dex. 



react strongly to M. All synthetic hydrogen Balmer line profiles 
show significant filling in due to wind emission for an increased 
mass loss, deteriorating the fit quality. Also the He i AAA! 1 line 
shows a decrease in core strength which is comparable to the 
decrease in the Ha core. This reconfirms that in order to self- 
consistently determine the mass loss rate all lines need to be 
fitted simultaneously. Therefore, a small discrepancy in the Ha 
core between the observed and synthetic line profile should not 
be considered a decisive reason to reject a fit. 

Exce pt fo r FHe the obtained parameters agree with the re- 
sults of lHPNl within the errors given by these authors. Similar 
to Cyg OB2 #8C we find a small helium enhancement. 

Cyg OB2 #4 The final fit to the spectrum of Cyg OB2 #4 is 
presented in Fig. [6] We obtained good fits for all lines, with 
exception of the helium singlet at 4922 A, for which the core is 
predicted too strong. However, recall that for this spectral type 



we assigned a relatively low weight to this line, for reasons 
explained in Sect. 13.21 

The parameters obtained from the fit agree well with the 
values of HPN, with exce ption of /3, for which we find a value 
higher by ~0.2. Note that HPN used a fixed value for fi to ob- 
tain their fit, whereas in this case the automated method self- 
consistently derived this parameter. 

The spectroscopic mass implied by the obtained logg 
value is significantly smaller than the evolutionary mass of 
Cyg OB2 #4. However, within the error bars (Sect. @J the two 
masses agree with each other. 

Cyg OB2 #1 In the final fit for this object, shown in Fig. 
there are two problematic lines. First, for the Hen /1468 6 line 
the core is predicted too strong. Even though compared to lHPNl 
the situation has improved considerably, the current version 
of fastwind still has difficulties predicting this line. Second, 
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Fig. 8. Same as Fig. [2] however for Cyg OB2 #2. Shown with dotted lines for Ha is the line profile of the best fit model with a M 
lower by a factor of 3. See text for further comments. 



the predicted Hei X441X is too weak. Possibly this is con- 
nected to the general ized dilution effect, for which we refer to 
iRepolust et alJJ200 4) for a recent discussion. 

In Fig. we also see that the Ha core of Cyg OB2 #10 
exhibits an emission feature. For this analysis we assumed that 
is was nebular and, consequently, excluded it from the fit. To 
test what the effect would be if this assumption was incorrect, a 
fit was made with this feature included in the profile. It turned 
out that the only parameter which was affected in this test was 
M, which showed a small increase of 0.04 dex. 



Cyg OB2 #2 For Cyg OB2 #2 the automated method could 
not self-consistently determine ft. Ther efore, we fixed its valu e 
at a theoretically predicted ft = 0.8 (cf. IPauldrach etalll986t> . 
In Fig. [S] the best fit is shown. We obtained good fits for all 
lines. However, in the case of Hei /14471 we do see a small 
under prediction of the forbidden component at 4469 A, which 
is likely related to incorrect line-broadening functions. 

For log g and M the obtained fit parameters differ consider- 
ably from the findings of HPN. We first focus on mass loss for 
which we obtain the relatively low rate of 1.63 x 10~ 7 M0yr _1 
with an error bar in the logarithm of this value of -0.15 and 
+0.12 dex (see Tab. given the quoted value of ft. Our M 
value is approxi mately a factor two higher than the mass loss 
rate obtained by HPN. These authors noted that it was not pos- 
sible to well constrain the mass loss rate of such a weak wind. 
Given the relatively modest errors indicated by our automated 
fitting method, we conclude that at least in principle our tech- 
nique allows to determine mass loss rates of winds as weak as 
that of Cyg OB2 #2. We have added the phrase "in principle" 
as it assumes the notion of beta and a very reliable continuum 
normalization, which is in this is case different from the one 
used by HPN for the He n 4541, He n 4686 and Ha lines. If this 
can not be assured, then systematic errors may dominate over 
the characteristic fitting error and the mass loss may be much 
less well constrained. Assuming the continuum location to be 



reliable, we illustrate the sensitivity of the spectrum to mass 
loss rates of ~10~ 7 M0yr~' by reducing the mass loss by a fac- 
tor of three. The Ha profile of this reduced mass loss model is 
shown in Fig.[8]as a dotted line. Comparison of these two cases 
shows that for winds of order 10~ 7 MQyr -1 the line still contain 
considerable M information. Interestingly, if we would not take 
into consideration the line core of Ha in our fitting method, we 
still recover the quoted mass loss to within three percent. We 
note that for our higher mass loss this object appears to behave 
well in the wind momentum luminosity relation (see Sect. 15. 41 . 
whereas IhPNI signal a discrepancy when using their estimated 
M value. 

The log g value obtained in this study is 0.36 dex larger than 
the value obtained "by eye' ' bv lHPNl Judging from the very 
good fits obtained, there is no indication that the automated fit 
overestimated the gravity. The spectroscopic mass of 17.0 Mq 
implied by the larger logg value, is also in good agreement 
with the evolutionary mass of Cyg OB2 #2, which is 18.7 Mq. 

3.3.2. Analysis of well studied dwarf OB-stars 

We have also reanalysed five well known and well studied 
dwarf OB-stars, sampling the range of O spectral sub-types, 
in order to probe a part of parameter space that is not well cov- 
ered by the Cyg OB2 stars. HD 217086 and £ Oph, for instance, 
are fast rotators with v r sin i = 350 and 400 km s _1 respectively 
(see also Tab.[5J- 10 Lac is a slow rotator, and r Sco is a very 
slow rotator. The latter two stars also feature very low mass loss 
rates, moreover, the actual M values of these stars are much de- 
H l2004 . HD 15629 is selected because 



bated (see Martin 
it appears relatively normal 



HD 15629 Apart from a slight over-prediction of the core 
strength in He n /14200, a very good fit was obtained for this 
object. The final fit is presented in Fig . |9] This object has re - 
cently been studied bv IRepolust et alJ J20041 hereafter IrPHIi . 
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Fig. 10. Same as Fig.|2] however for HD 217086. 



Compared to the parameters obtained by these authors, we find 
good agreement except for T e g, M and p. Note that in contrast 
to this study we do not find a helium deficiency. However, the 
difference of 0.02 with respect to the solar value obtained here 
is within the error quoted by RPH. 

The difference in wind par amete rs can be explained by the 
value of (3 — 0.8 assumed by IrphL Our self-consistently de- 
rived value for (3 — 1.18. As the effect of (3 on the spectrum is 
connected to the mass loss rate through the velocity law and the 
continuity equation, the lower M obtained with the automated 
method is explained. 

The 1 .5 kK increase of compared to [RPH! can be at- 
tributed to the improved fit quality and the increase in log g of 
0.1 dex. An increase in logg implies an increase in electron 
density, resulting in an increase in the recombination rate. The 
strength of both the He i and He n lines depend on this rate, as 
the involved levels are mainly populated through recombina- 
tion. Consequently, as He m is the dominant ionization stage in 



the atmosphere of HD 15629 the strength of the He i and He n 
lines will increase when the recombination rate increases. To 
compensate for this increase in line strength an increase in T e g, 
decreasing the ionization fractions of He i and He n, is neces- 
sary. 

HD 21 7086 With a projected rotational velocity of 350 km s _1 
this object can be considered to be a fast rotator, and our analy- 
sis of this object will show how well the automated method can 
handle large v r siru. In Fig. [TO] the best fit obtained with our 
method is presented. We find that the large projected rotational 
velocity does not pose any problem for the method, i.e. the fit 
quality of all the lines fitted is very good. 

With respect to the obtained parameters, again, these can 
be compared to the work of RPH. In this comparison we find 
considerable differences for T e g and \ogg and a small differ- 
ence for Tee - The effective temperature found by the automated 
method is 2.1 kK higher. This is a significant increase, but 
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when the log g value obtained here is considered, this can be 
explained in a similar manner as the T s g increase of HD 15629. 

The best fit is obtaine d with a log g value that is 0.29 dex 
higher than the value from RPH. Judging from the line profiles 
in Fig. 1 101 there is no evidence for an overestimation of logg. 
T his higher log g remov es the disc repanc y with the calibr atio n 
of iMarkova et all (2004) found bv lRPHJ (see Fig. 17 in lRPHl) . 
We also note that, similar to Cyg OB2 #2, the increased log g 
implies a spectroscopic mass which agrees well with the evo- 
lutionary mass of HD 217086 (cf. Tab.0}. This is not the case 
for the value determined bv lRPHl which points to a clear dis- 
crepancy. 

The considerable helium abundance enhancement found by 
IRPHI is not reproduced by the automated method. Even though 
this object is a rapid rotator, our fit indicates a normal, i.e. solar, 
helium abundance. 



1 Lac Like in the case of Cyg OB2 #2, and for the remaining 
objects, the wind is too weak to self-consistently determine (3. 
Therefore, again a value of/3 - 0.8 was assumed. 

The photospheric parameters obtained for 10 Lac agree 
very well with the results of IHPNI The b est fit to the ob- 
served spectrum is shown in Fig. ^2 Whereas HPN find that the 
mass loss rate cannot be constrained and only an upper limit of 
10~ 8 M0yr _1 is found, the automated method was able to self- 
consistently determine M at 6x 10~ 8 Mgyr -1 , though with large 
error bars (see Tab.|5}. Our error bar indicates that M may be 
an order of magnitude lower, i.e. it may still be consistent with 
the iHPNl result. 

Various other authors have determined the mass loss rate 
of 10 Lac using different method s. These determ inations range 
from up to 2 x 10~ 7 f Howarth & Prinia 1989) down to 2 x 
10~ 9 M yr-' (Mar tins et al.ll2004 . Consequently, compared 
to these independent determinations no conclusive answer can 
be given to the question whether the M derived from the opti- 
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cal spectrum is correct. We conclude that the mass loss rate of 
10 Lac is anomalously low when placed into context with the 
other dwarfs stars studied here. For instance, the dwarfs f Oph 
and HD 217086, which have luminosities that are, respectively, 
lower and higher by ~0.1 dex, both exhibit a mass loss rate 
higher by several factors. In Sect. I5.4l we will discuss this fur- 
ther in terms of the wind-momentum luminosity relation. 



£ Oph The large v r sku of 400 kms" 1 was not a problem to 
obtain a good fit. In Fig. El the best fit for f Oph is pre- 
sented. With the exception of the helium abundance, the com- 
parison with the results of RPH yields very good agreement. 
Note that the mass loss rate obtained by these authors is an 
upper limit, whereas in this study M could be derived self- 
consistently. With respect to iHe we do not find any evidence 
for a significant overa bundance of helium, in agreement with 
IViUamariz & Herrerdd2005h . 



t Sco The best fit for t Sco is presented in Fig.^] All lines, 
including H6, which is not shown here, are reproduced accu- 
rately. The photosp heric parameters we obtai ned can be com- 
pared to the work of lSchonberner et all dl988l) and lKilian et alJ 
(119911) who b oth studied t Sco using plane parallel mod- 
els. iKilian et alJ fo und r e ff=31.7 kK and log,g=4.25, whereas 
ISchonberner et al.l obtained r e ff=33.0 kK and logg=4.15. The 
difference in T e ff between the two studies is explained by the 
fact that in the latter analysis no line blanketing was included 
in the models. Therefore, we prefer to compare our T s g to the 
former investigation, which agree very well. In terms of the 
gravity we find good agreement with the second study. The 
value obtained bv lKilian et alJ seems rather high. Given the al- 
most perfect agreement between the synthetic line profiles and 
the observations in Fig. [D] the reason for this discrepancy is 
unclear. On a side note, more recently Repol ust et a 1. (2005) 
analysed the infrared spectrum of this object. Their findings do 
confirm our lower value, but could not reproduce the enhanced 



helium abundance we find, due to a lack of observed infrared 
He ii lines. 

In the recent literature the mass loss rate usually adopted 
for t Sco is 9 x 10~ 9 Mgyr -1 , which is considerably smaller 
than the 6.14 x 10~ 8 Moyr -1 obtained in this study. However, 
the former mass loss rate is an average value determined by 
Ide Jager et alJ (1988), based on the mass loss rates indepen- 
dently found by Gathi er et alJ (11981) and Hamann dl98lli . 
Based on the UV resonance lines, these two studies, respec- 
tively, determined M to be 7.4 x 10~ 8 and 1.3 x 10~ 9 M Q yr'. 
So, they differ by more than a factor of 50. The mass loss rate 
obtained with the automated method is in reasonable agreement 
with that obtained by Gathi er et all Our higher value is also 
supp orted by the stu dy of the infrared spectrum of t Sco by 
Rep olust etai1ll2005l) whofindM 2xl0~ 8 M Q yr _1 . Detailed 
fitting of Bra will likely clarify this issue. 



4. Error analysis 

Here we will introduce our method of estimating errors on the 
parameters derived with the automated method. This method 
is based on properties of the distribution of the fitnesses of the 
models in parameter space, which may seem conceptually dif- 
ferent from classical approaches of defining error bars (and in 
a sense it is). However, we will demonstrate for the case of 
10 Lac that our error definition is very comparable to what is 
routinely done in fit diagram approaches. 

4.1. Fit diagrams 

In a fit diagram method the error bar on r e ff and log g is de- 
rived by investigating the simultaneous behaviour of these two 
parameters. In panel a of Fig.[H]the fit diagram of 10 Lac is 
presented adopting for all other parameters (save T e g and log g) 
the best fit values obtained in our automated fitting. This dia- 
gram was constructed by calculating a grid of fastwind models 
in the r e ff-log g plane, and evaluating for every line for every 
r e ff which model, i.e. logg, fits this line the best. The loca- 
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Fig. 14. Panel a: fit diagram of T e ff and \ogg for 10 Lac. Panel b: Fitness as a function of for logg = 4.0. Panel c: Fitness 
distribution of the models calculated during the fitting run of the automated method. Panel d: Distribution of in the models 
located within the global optimum. The maximum variation of within the global optimum, which corresponds to the error 
estimate of this parameter, is ~900 K. 



tion where the resulting fit curves intersect, corresponds to the 
best fit. This best fit yields r eff = 36000 K and logg = 4.0. 
Note that this result was obtained without the use of our auto- 
mated method. The error can now be estimated by estimating 
the dispersion of the fit curves around this location. In panel a 
of Fig. ^]this is indicated by a box around the best fit loca- 
tion. The corresponding error estimates are 1000 K in and 
0.1 dex in logg. 



The method described above cannot be applied to our au- 
tomated fitting method due to two reasons. First, as we have 
defined the fit quality according to Eq. ([0, this definition of 
fitness compresses the fit curves of all individual lines in the 
fit diagram to a single curve. In Fig.[H]this curve is shown as 
a thick dashed line. Although the curve runs through the best 
fit point, no information about the dispersion of the solutions 
around this point can be derived from it. The second reason 
lies in the multidimensional character of the problem of line 
fitting. If one would want to properly estimate the error tak- 
ing this multidimensionality in to account, a fit diagram should 
be constructed with a dimension equal to the number of free 
parameters evaluated. In case of our fits this translates to the 
construction of a six dimensional fit diagram. 



4.2. Optimum width based error estimates 

Even though we have argued that fit diagrams cannot be used 
with our fitting method, it is possible to construct an error esti- 
mate which is analogous to the use of these diagrams and does 
take the multidimensionality of the problem into account. This 
can be done by first realizing that the error box shown in Fig. 1141 
essentially is a measure of the width of the optimum in parame- 
ter space, i.e. it defines the region in which models are located 
which approximately have the same fit quality. This is illus- 
trated in panel b. There we show the one-dimensional fitness 
function in the r e ff-log g plane for log g = 4.0. Indicated with 
dashed lines is the error in r e ff estimated using the fit diagram 
of 10 Lac. Confined between these lines is the region which 
corresponds to the optimum as defined by the error box in panel 
a. Consequently, the difference between maximum and mini- 
mum fitness in this region defines the width of the optimum. 
Returning to the general case, we can now invert the reason- 
ing and state that the error estimate for a given parameter is 
equal to the maximum variation of this parameter in the group 
of best fitting models, i.e. the models located within the error 
box. Consequently, in the automated fitting method by defining 
a group of best fitting models, the error estimates for all free 
parameters can be determined. 

We define the group of best fitting models as the group 
of models that lie within the global optimum. Put differently, 
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Table 5. Error estimates for fit parameters obtained using the automated fitting method and parameters derived from these. 
Denoted by nd are errors in v tur b that reach up to the maximum allowed value of v tur b and, therefore, are formally not defined. 
Uncertainties in the fit parameters result from the optimum width based error estimates method. See text for details and discussion. 
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the width of the global optimum in terms of fitness, defines 
the group of best fitting models. Identifying and, consequently, 
measuring this width is facilitated by the nature of the GA, i.e. 
selected reproduction, incorporated in our fitting method. Due 
to this selected reproduction the exploration through parame- 
ter space results in a mapping of this space in which regions 
of high fit quality, i.e. the regions around local optima and the 
global optimum, are sampled more intensively. Consequently, 
if we would rank all models of all generations calculated dur- 
ing a fitting run according to their fitness, the resulting distri- 
bution will peak around the locations of the optima. In case of 
the global optimum the width of this peak, starting from up to 
the maximum fitness found, is, analogous to the width of the 
error box used in a fit diagram, a direct measure of the width of 
the optimum. Consequently, this width depends on the quality 
of the data, i.e. it will be broader or narrower for, respectively, 
low and high signal to noise, and on the degeneracy between 
the fit parameters. Therefore, the error estimates of the individ- 
ual parameters are equal to the maximum variations of these 
parameters for all models contained in the peak corresponding 
to the global optimum. 



In panel c of Fig.[H]the distribution of the models accord- 
ing to their fitness calculated during the fitting run of 10 Lac 
using the automated method is shown. The fitnesses are nor- 
malized with respect to the highest fitness and only the top half 
of the distribution is shown. In this distribution two peaks are 
clearly distinguishable. The most pronounced peak is located 
at F w 0.9 and corresponds to the region around the global op- 
timum. A second peak, corresponding to a region around a sec- 
ondary optimum, is located at F « 0.83. To derive the error on 
the fit parameters we estimate the total width of the global op- 



timum for 10 Lac to be -0.15 2 , i.e. the range of F — 0.85. ..1.0 
corresponds to the width of the optimum. In panel d of Fig. 1141 
we show the resulting distribution of T e ff of the models within 
this global optimum. In this figure we see that the maximum 
variation, hence the error estimate, is -900 K, which is in 
good agreement with the value derived using the fit diagram 
of 10 Lac. For \ogg we also find an error estimate of -0.1 dex, 
which is also very similar to the value obtained with this di- 
agram. The exact values as well as error estimates for all fit 
parameters of all objects are given in Tab. [5] It is important 
to note that our error analysis method also allows for an error 
estimate of parameters to which the spectrum does not react 
strongly. For 10 Lac this is clearly the case for the mass loss 
rate, for which we find large error bars. 

4.3. Derived parameters 

In Tab.[5]the errors on the derived parameters were calculated 
based on the error estimates of the fit parameters. Here we will 
elaborate on their derivation. 

The error in the stellar radii is dominated by the un- 
certainty in the absolute visual magnitude. In case of the 
Cyg OB2 objects we adopt these to be 0.1 m conform the work 
of Massev & Thompsonl Jl99lb . For HD 1562 9, HP 217086 
and f Oph we use the uncertainty as given by RPHJ of 0.3 m . 
The distance to 10 Lac and r Sco was measured by Hipparcos. 
Therefore, for these two objects we adopt the error based 
on this measurement, which, respectively, is 0.4 and 0.2 m . 
Together with the uncertainty in r e +> the uncertainty in R+ is 



2 In general this is not a fixed number. Considering all programme 
stars we find the width of the global optimum to be within the range 
-0.1 to -0.2. 
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Fig. 15. Comparison of the effective temperatures obtained us- 
ing automated fits (horizontal axis) and "by eye" fits. On the 
vertical axis the ratio of automated relative to "by eye" tem- 
perature determination is given. The dashed lines correspond 
to a four percent error usually adopted for "by eye" determined 
values. 




Fig. 16. Gravities obtained with automated fits (horizontal axis) 
are compared to gravities determined from "by eye" fits. The 
vertical axis gives the difference of the logarithm of the two 
gravity determinations. Indicated by dashed lines are the 0.1 
dex errors usually adopted for gravities determined "by eye". 



calculated according to Eq. (8) of RPH, where we used the 
largest absolute uncertainty in T e ff for a given object. 

To correct the surf ace gravity for centrifugal forces, a cor- 
rection conform Herrer o et al.l I I 19921) was applied to the gravity 
determined from the spectral fits. This corrected value is given 
in Tab. 0] As shown by RPH this correction has a non negligi- 
ble effect on the error in the resulting log g c . Consequently, we 
used their estimate to calculate the total error estimate of log g c 
as given in Tab. [5] Using this error together with the uncertainty 
in R+ the resulting uncertainty in the spectroscopic mass was 
calculated. 

For the calculation of the uncertainty in the stellar luminos- 
ity, we consistently adopted the largest absolute error in r e g\ 
The resulting A log L* as well as the uncertainty in T e g have 
an effect on the evolutionary mass. We have estimated errors 
for this quantity using the error box spanned by A log L* and 
A log T eS . 



5. Comparison with previous results 

In this section we will compare the results obtained with our au- 
tomated fitting method with those from "by eye" fits (relevant 
references to the comparison studies are given in the previous 
section). This does not constitute a one-to-one comparison of 
the automated and "by eye" approach as this would require the 
use of identical model atmosphere codes as well as the same set 
of spectra, moreover, with identical continuum normalization. 
Potential differences can therefore not exclusively be attributed 
to the less bias sensitive automated fitting method. However, as 
we have applied our method to a sizeable sample of early type 
stars, the automated nature of it does assure that it is the most 
homogeneous study to date, i.e. without at least some of the 
biases involved in conventional analyses. 



5.1. Effective temperature 

In Fig. [B] a comparison of the effective temperatures deter- 
mined in this study with r e ff values obtained with "by eye" fits, 
is presented. Indicated with dashed lines are the four percent 
errors usually adopted for "by eye" fitted spectra. With the ex- 
ception of the outlier HD 217086 at 38.1 kK, the agreement is 
very good and no systematic trend is visible. From this plot we 
can conclude that the T e g obtained with the automated fit is at 
least as reliable as the temperatures determined in the conven- 
tional way. 

5.2. Gravities 

In many cases the gravity obtained with the automated proce- 
dure is significantly higher than the values obtained with the 
conventional "by eye" fitted spectra. This is shown in Fig. 1161 
where we show as a function of the gravities obtained in this 
study the differences with the "by eye" determined values. 
Indicated with dashed lines in this figure is the 0.1 dex error 
in logg that is often assigned to a "by eye" fitting of the hy- 
drogen Balmer line wings. It is important to note that this plot 
shows that there is no obvious trend in the differences, i.e. there 
appears no systematic increase as a function of log g present. 

It is clear, however, that there are three outliers for which 
previous gravity determinations yield values that are at least 
0.2 dex lower. These are in order of increasing gravity (as 
determined in this study): Cyg OB2 #2, Cyg OB2 #7 and 
HD 217086. For all three cases, previous spectroscopic mass 
determinations result in values that are about a factor of two 
less than the corresponding evolutionary masses. One reason 
for these discrepant gravity values can be traced to a difference 
between automated and "by eye" fitting. In "by eye" fitting, it 
is custom to prohibit the theoretical line flux in the wings of 
Balmer lines - specifically that at the position in the observed 
line wing where the profile curvature is maximal - to be be- 
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low th at of the observed flux. This constraint has been used by 
iHPNl for the Cyg OB2 stars; for HD 217086, we could not ver- 
ify whether this was the case. The automated method does not 
apply this constraint. Therefore, as it strives for a maximum fit- 
ness, it tends to fit the curve through the signal noise as much 
as possible. This yields a higher gravity. 

A second reason is connected to the multidimensional na- 
ture of the optimization problem. "By eye" fitting may not find 
the optimum fit, as in general it can not simultaneously deal 
in a sufficiently adequate way with all the free parameters of 
the problem. Consequently, some of the "by eye" fitted spec- 
tra do not correspond to the best fit possible. A good example 
in which this appears to be the case is HD 217086. With the 
automated fit we not only obtained a gravity that is higher by 
~0.3 dex, but also an effective temperature higher by 2.1 kK 
compared to the results o f RPH. Consequently, as the ioniza- 
tion structure of the atmosphere depends heavily on this tem- 
perature, so does the gr avity one obtains from a spectral fit 
for this temperature. As RPH obtained a gravity for a signif- 
icantly lower effective temperature, the gravity obtained from 
their spectral fit likely corresponds to the value from a local 
optimum in parameter space. 

5.3. Helium abundance and microturbulence 

This analysis is the first in which the helium abundance and the 
microturbulent velocity have been tre ated a s continuous free 
parameters. In the studies of iHPNl and lRPHl only two possible 
values for the microturbulent velocity were adopted. For the he- 
lium abundance an initial solar abundance was adopted, which 
was modified when no satisfying fit could be obtained for this 
abundance. Consequently, a comparison with these studies as 
was done for e.g. the gravities, is not possible. Instead we will 
only discuss whether the obtained values of these parameters 
are reasonable and comment on possible correlations with other 
parameters. 

The helium abundances given in Tab. 0] show that no ex- 
treme values were needed by the fitting method to obtain a 
good fit. An exception to this may be Fee =0.21 obtained for 
Cyg OB2 #7. However, as discussed earlier this value i s still 
significantly smaller than the iHe=0.3 obtained bv lHPNl With 
respect to a possible relation between the helium abundance 
and other parameters, only a small correlation between T e s and 
Tne is found for the supergiants. For these objects it appears 
(cf. Tab. 0} that the helium abundance increases with increas- 
ing effective temperature. However, as we only analysed six 
supergiants further investigation using a larger sample needs to 
be undertaken. 

Also for the microturbulent velocities no anomalous values 
were needed to fit the spectra. The large error bars in the turbu- 
lent velocity quoted in Tab. [5] especially for the supergiants, 
show that the profiles are not very sen sitive to this parame- 
ter. Th is is consistent with the study of Villamariz & Herrero 
(2000) and lRPHl The nd entries given on some of the positive 
errors in the table indicate that they reach up to the maximum 
allowed value of v tur b, which is 20 kms -1 . Therefore, they are 
formally not defined. The fact that some of the small scale tur- 
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Fig. 17. Difference between mass loss rate obtained by the au- 
tomated method (given by the horizontal) axis and values de- 
termined by eye. A typical 0.15 dex error is indicated by the 
dashed lines. The two outliers at log M -7.2 and and log M 
-6.8, respectively, correspond to 10 Lac and Cyg OB2 #2. 



bulent velocities are close to this maximum value may indicate 
that they represent lower limits, though, again, this likely re- 
flects that they are poorly constraint. 

No correlation of the microturbulence with any of the other 
parameters is found. In particular not between v tur b and log g 
and v tur b a nd F Hp . Vario us authors have hinted at such a corre- 
lation (e.g. lKilianll992l 

5.4. Wind parameters 

The straightforward comparison of the mass loss rates obtained 
with the automated method with values determined from spec- 
tral fits "by eye" is shown in Fig. El With exception of t Sco 
at logM = -7.2, f or which the mass loss rate determined by 
iGathier et alJ J 1 98 ll) from UV line fitting serves as a compar- 
ison, all mass loss rates are compared to values determined 
from Ha fitting. For this comparison we assume an error of 
0.15 dex in the "by eye" determined values. This uncertainty 
corresponds to a typical error obtained from Ha fitting and is 
shown in Fig.^l as a set of dashed lines. With exception of 
10 Lac and r Sco, for which the mass loss rate determination 
is uncertain, this error is also comparable to the errors obtained 
with the automated method. 

Two objects show a relative increase in M which is much 
larger than the typical error. These are 10 Lac at log M -7.3 
and Cyg OB2 #2 at log M -6.8. In the case of the latter we 
showed that the increase is due to a more efficient use of wind 
information stored in the line profiles by the automated method, 
which improves the relation of Cyg OB2 #2 with respect to the 
wind-momentum relation (see Sect. 16. 21 . 

With respect to 10 Lac we already mentioned that a range 
of more than two orders of magnitude in mass loss rate has 
been found in different studies. Her e we h ave made the com- 
parison with the upper limit found bv lHPNl which corresponds 
to one of the lowest M determined for this object. If we would 
have compared our findings to the higher value obtained by 
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Fig. 18. Spectroscopic masses derived in thi s study compared 
to evolutionary masses from Schalle r et alJ II1992I) . With the 
gravities obtained from the automated fits no mass discrepancy 
is found and no systematic deviation between the spectroscopi- 
cally derived masses and the evolutionary predicted masses can 
be observed. 
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Fig. 19. Modified wind momentum (MWM) in units of 
[gems -2 /?©] °f me objects fitted with the automated method 
(solid dots). The solid line, giving the wind-momentum lumi- 
nosity relation (WLR), corresponds to the regression of the 
modified wind momenta. Given by the dashed line is the pre- 
dicted WLR of lVink et alJ i2000fc . 



iHowarth & Prinial Jl989h . 10 Lac would be at AM = -0.5, i.e. 
the situation in Fig.Elwould be reversed. Consequently, the 
large difference for 10 Lac shown in this figure can not be as- 
signed to an error in the automated method, but rather reflects 
our limited understanding of this object. 

All in all, we can conclude that the general agreement be- 
tween mass loss rates obtained with the automated method and 
"by eye" determinations is very good. 

6. Implications for the properties of massive stars 

With our automated method we have analysed a sizeable sam- 
ple of early type stars in a homogeneous way, which allows 
a first discussion of the implications the newly obtained pa- 
rameters may have on the mass and modified wind-momentum 
luminosity relation (WLR) of massive stars. A thorough dis- 
cussion however needs to be based on a much larger sample, 
therefore at this point we keep the discussion general and the 
conclusions tentative. 

6.1. On the mass discrepancy 

The so call ed mass discrepancy was first noticed by 

iHerrero et alJ Jl992h . These authors found that the spectro- 
scopic masses, i.e. masses calculated from the spectroscopi- 
cally determined gravity, were systematically smaller than the 
masses predicted by evolutionary calculations. The situation 
improved considera bly with the use o f unified stellar atmo- 
sphere rrwdd^J^gjHe^^^^OllOQ^)- However, as pointed 
out by Repol ust et alJ (I2004h for stars with masses lower than 
50 Mq still a milder form of a mass discrepancy appears to 
persist. 

Does the automated fitting method, employing the latest 
version of fastwind, help in resolving the mass discrepancy? 
In Fig.[2Dwe present a comparison of the spectroscopic masses 



calculated with the gravities obtained in this study, with masses 
derive d by interpolating evolutionary tracks of ISchaller et alJ 
d!992h . It is clear that with the new gravities the situation is 
very satisfying. All objects have spectroscopic and evolution- 
ary masses which agree within the error bars. 

For stars with masses below 50 Mq a milder form of the 
mass discrepancy (as found bv iRPHt see their Fig. 20) could 
still be present, but with the present data no systematic offset 
between the two mass scales can be appreciated. Though we 
feel it may be premature to conclude that the present analysis 
shows that the mass discrepancy has been resolved, our results 
point to a clear improvement. 

6.2. Wind-momentum luminosity relation 

The modified stellar wind momentum (MWM) versus luminos- 
ity relation offers a meaningful way to compare observed wind 
properties with as pects and predictions of the theory of line 
driven winds (see iKudritzki & PulslEoO O for a comprehensive 
discussion). Without going into any detail, the modified wind 

1 II 

momentum D mom - Mv x R^ is predicted to be a power law of 
stellar luminosity. 

logD mom = xlog(L±/L©) + log£>o , (7) 

where x, the inverse of the slope of the line-strength distribu- 
tion function corrected for ionization effects (IPuls et al. 2000), 
is expected to be a function of spectral type and metal abun- 
dance, and Dp is a function of m etallicity and possibly lumi- 
nosity class dMarkova et al.l2004l) . The advantageous property 
of D mom is that it is not very sensitive to the stellar mass. 

The limited number of stars studied in this paper is clearly 
insufficient to disentangle subtleties in the D mom vs. L+ rela- 
tion. However, it is interesting to compare the observed and 
predicted modified wind momentum, as well as to discuss the 
location of 10 Lac - a notorious outlier. 
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Figure [23 shows this comparison between derived and the- 
oretical modified wind momentum. Using all programme stars 
to construct an empirical linear curve in the units of this dia- 
gram gives the following relation 

logD mom = (1.88 ± 0.09)log(L*/io) + (18.59 + 0.52) . (8) 

Within the given erro rs this relation is equal to the theoreti- 
cal WLR predicted bv lVink et alJ J2000h . who found x = 1.83 
and logD mom = 18.68. Note that the low luminosity objects 
(log L^/Lq < 5.5) also follow the average relation. Therefore, 
the newly obtained mass los s rat es do not show the discrep ancy 
found bv lPuls et all d 1996b and iKudritzki & Pulsl (|2000), but 
confirm the work of RPH. These authors found the low lumi- 
nosity objects to follow the general trend, based on upper limits 
they obtained for the mass loss rates, whereas our new method 
is sensitive enough to determine these self-consistently. 

In order to investigate the effect of the anomalously low 
D mom obtained for 10 Lac, we also constructed a WLR exclud- 
ing this object. We found that for this new relation the parame- 
ters x and log D only changed with ~0.02 and ~0.01, respec- 
tively, reflecting the large erro r bars found for this ob ject. 

Previous investigations by Mark ova et alJ d2004h and lRPlj 
have found the WLR to be as function of luminosity class. 
Whereas the former study finds a steeper WLR for the su- 
pergiants compared to the dwarfs, the latter finds the opposite 
(though iRPHl remark that the subset of Cyg OB2 stars seem 
to behave more in accordance with the theoretical result). In 
our sample no obvious separation is visible. In particular note 
the two objects overlapping at log = 4.9 in Fig. [HJ] 
which are the dwarf £ Oph and the supergiant Cyg OB2 #2. 
To investigate a possible separation in more detail, a separate 
WLR was constructed for the Cyg OB2 supergiants. The re- 
sulting values of the parameters obtained are x = 1.79 + 0.14 
and log D = 19.12 + 0.80. The d ecrease in x qualitatively con- 
firms the work of Mar kova et all However, we have to realize 
that our sample might be too small from a statistical point of 
view to be able to firmly conclude whether a real separation 
exists. Therefore, this question has to be postponed until we 
have analysed a larger sample. 

7. Summary, conclusions and future work 

We have presented the first method for the automated fitting of 
spectra of massive stars with stellar winds. In this first imple- 
mentation, a set of continuum normalized optical spectral lines 
is fitted to predictions made with the fast perform ance non-LTE 
model atmosphere code fastwind bv lPuls et alJ d2005l) . The fit- 
ting method i tself is based on the genetic algorithm pikaia by 
Charbonnei ^] Jl995l) . which was parallelized in order to handle 
the thousands of fastwind models which have to be calculated 
for an automated fit. Concerning the automated method we can 
draw the following conclusions: 

i) The method is robust. In applying the method to a num- 
ber of formal tests, to the study of seven O-type stars in 
Cyg OB2, and to five Galactic stars including extreme rota- 
tors and/or stars with weak winds (few times 10~ 8 M0yr~') 
the fitting procedure did not encounter convergence prob- 
lems. 



ii) Using the width of the global optimum in terms of fitness, 
defining the group of best fitting models, we are able to de- 
fine error estimates for all of the six free parameters of the 
model (T e ff, logg, helium over hydrogen abundance, v tur b, 
M and /?). These errors compare well with errors adopted 
in "by eye" fitting methods. 

Hi) For the investigated dataset our automated fitting method 
recovers mass-loss rates down to ~ 6 x lO~ 8 M0yr~' to 
within an error of a factor of two. We point out that even 
for such low mass-loss rates it is not only the core of the 
hydrogen Ha line that is a mass-loss diagnostics. When ig- 
noring this core the GA still recovers M, showing that the 
GA is also sensitive to indirect effects of a change in M 
on the atmospheric structure as a whole. However, for the 
method to fully take advantage of this information a very 
accurate continuum normalization is required. 

iv) Though we have so far tested our method for O-type stars 
and early B-type dwarf stars, the method can also be ap- 
plied to B and A supergiants when atomic models of diag- 
nostic lines (such as Si in and Si iv) are implemented into 
the analysis. 

We have re-investigated seven O-type stars in the young clus- 
ter Cyg OB2 and compared o ur re sults with the study by 
iHerrero. Puis. & Naiarrol d2002h . The iHPNl study uses an ear- 
lier version of fastwind and a "by eye" fitting procedure. The 
only difference between the two studies in terms of the treat- 
ment of the free parameters is that HPN did not treat the mi- 
croturbulent velocity and the hydrogen over helium abundance 
ratio as continuous free parameters. Instead, they opted to adopt 
in case of the former two possible values only. In case of the 
latter an initial solar abundance was adopted which was mod- 
ified in case no satisfying fit could be obtained with this solar 
abundance. We have also compared the results of an automated 
fitting to five early-type dwarf stars to further investigate the 
robustness of our method for stars with high rotational veloci- 
ties and/or low mass loss rates. With respect to weak winds we 
refer to conclusion Hi. Regarding large v r sin i values, we find 
that these do not pose problems for the automated method. This 
is reflected in conclusion i. Concerning the spectral analysis of 
the entire sample we can draw the following conclusions: 

v) For almost all parameter s we fi nd excellent agreement with 
the results of IHPNl and IRPHl which, we note, make use 
of a previous version of fastwind and an independent con- 
tinuum normalization. The quality of our fits (in terms of 
fitness, which is a measure for the \ 2 of the lines) is even 
better than obtained in these prior studies. 

vi) In three cases we find a significantly higher surface gravity 
(by up to 0.36 dex). We identify two possible causes for this 
difference that may be connected to the difference between 
automated and "by eye" fitting. First, comparison of the two 
methods indicates that in fitting the Balmer line wings the 
latter method places essentially infinite strength to the ob- 
served flux at the point of maximum curvature of the wing 
profile. The automated method does not do this. Second, as 
the automated method is a multidimensional optimization 
method it may simply find a better fit to the overall spec- 
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trum. In at least one case this implied a higher temperature 
and significantly higher gravity. 

vii) A comparison of our derived masses with those predicted 
by evolutionary calculations does not show any system- 
atic discrepancy. Su ch a discrepancy was first noted by 
Herrero et alJ Jl992l) . though was part ly res olved when 
model atmospheres improved (e.g. see 1HPN1) . Still, with 
state-of-the-art models a mild form of a mass discrepancy 
remained for stars with masses below 50 Mq (e.g. see 
iRPHh . The automated fitting approach in combination with 
the improved version of fastwind does not find evidence for 
a mass discrepancy, although we remark that a truly robust 
conclusion, particularly for stars between 20 and 50 solar 
masses, may require the investigation of a larger sample. 

viii) The empirical modified wind momentum relation con- 
structed on the basis of the twelve objects analysed in this 
study agree to within the e rror bars with the theoretical 
MWM relations based on the lVinketal] 12000 ) predictions 
of mass loss rates. 

This first implementation of a genetic algorithm combined 
with the fast performance code fastwind already shows the 
high potential of automatic spectral analysis. With the current 
rapid increase in observations of early-type massive stars the 
need for an automated fitting method is evident. We will first 
use our method to analyse the ~100 O-type and early B-type 
stars observed in t he VLT large prog ramme FLAMES Survey 
of Massive Stars ( Eva ns et alJ 120051) in the Galaxy and the 
Magellanic Clouds in a homogeneous way. Future develop- 
ment of the automated fitting method is likely to be in conjunc- 
tion with the further development of fastwind. Improvements 
will include the mo deli ng of: near-infrared lines (see e.g. 
iLenorzer et alJl2004l and iReoolust et alJ Eo05). optical CNO 
lines fsee e.g. lTrundle et all2004l) . and possibly UV resonance 
lines. Additional model parameters that may be constrained 
within an automated approach include a depth dependent pro- 
file for the microturbulent velocity and small scale clumping. 
Within the current implementation, most likely the method will 
also be able to constrain the terminal flow velocity of A -type 
supergiants jMcCarthv et al.ll997tlKudritzki et alll999t) . 
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